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Abstract-  As  the  Navy  expands  its  dependence  on  underwater 
communication  between  sensors,  operating  in  areas  where 
turbulent  turbid  water  layers  are  present;  there  exists  a  need  to 
accurately  predict  prior  to  sensor  deployment  how  they  will 
operate  in  these  environments.  The  benthic  nepheloid  layer  (BNL) 
is  an  example  of  a  moving  turbid  water  layer  in  the  ocean.  The 
BNL  is  characterized  by  changing  vertical  thickness, 
concentration  and  speed  of  the  suspended  material.  Acoustic 
propagation  and  hence  acoustic  communication  and  system 
performance  will  be  affected  when  operating  in  these  areas.  The 
suspended  material  will  alter  the  sound  speed,  density  and  the 
attenuation  of  the  medium.  Thus  what  was  once  a  non-dispersive 
quiescent  environment  is  now  a  moving  dispersive  environment. 
The  numerical  solution  of  acoustic  pulse  propagation  through 
dispersive  moving  media  requires  the  inclusion  of  attenuation 
and  its  causal  companion,  phase  velocity.  For  acoustic 
propagation  in  a  linear  dispersive  quiescent  medium,  Szabo  [J. 
Acoust.  Soc.  Am.,  96,  491-500  (1994)]  introduced  the  concept  of  a 
convolutional  propagation  operator  that  plays  the  role  of  a  casual 
propagation  factor  in  the  time  domain.  The  operator  has  been 
incorporated  in  the  linear  wave  equation  for  quiescent  media. 
Additionally  it  has  been  used  to  study  propagation  and  scattering 
from  such  widely  diverse  media  as  bubble  plumes  in  the  ocean 
and  ultrasound  propagation  in  human  tissue.  In  this  work,  this 
method  is  extended  to  address  acoustic  propagation  in  dispersive 
moving  media.  The  development  of  the  modified  wave  equation 
for  sound  propagation  in  dispersive  media  with  inhomogeneous 
flow  will  be  described,  along  with  an  example.  The  resulting 
modified  wave  equation  is  solved  via  the  method  of  finite 
differences. 

I.  Introduction 

Traditionally  underwater  acoustic  propagation  assumes  a 
quiescent  environment.  That  is,  since  the  acoustic  signal  in 
sea- water  travels  nominally  at  1500  m/s,  the  travel  time 
between  the  source  and  some  arbitrary  receiver  is  less  than,  or 
is  assumed  to  be  less  than,  the  time  it  takes  the  environment  to 
change  appreciably  enough  such  that  the  acoustic  field 
experiences  this  change.  There  are  however  situations  where 
this  assumption  doesn't  hold.  For  example  when  a  benthic 
nepheloid  layer  (BNL)  is  present.  A  BNL  is  formed  when 
sediments  lying  on  the  ocean  floor  are  resuspended  due  to 
water  motion.  The  water  could  be  in  motion  due  to  currents, 
tides  or  via  the  shoaling  and  breaking  of  internal  waves  on  the 
continental  shelf.  The  BNL  is  dynamic,  characterized  by 
changing  vertical  thickness,  concentration  and  speed  of 
resuspended  material.  This  increase  in  concentration  of 
material  suspended  in  the  water  column  creates  a  dispersive 


environment.  When  the  environment/medium  is  dispersive, 
causality  dictates  that  the  propagating  field  suffers  from 
attenuation.  The  attenuation  is  frequency  dependent.  In 
addition  to  this  frequency  dependent  attenuation  is  the  causal 
phase  velocity.  When  the  acoustic  source  is  time  limited,  i.e. 
an  acoustic  pulse,  it  is  composed  of  many  frequencies.  Thus  to 
accurately  model  such  a  situation  requires  that  the  solution  be 
able  to  take  into  account  the  appropriate  attenuation  and  phase 
velocity  of  all  frequencies  making  up  the  source  signature. 
Direct  modeling  in  the  time  domain  is  highly  desirable,  since 
it  allows  for  more  direct  and  efficient  numerical  solutions,  and 
causality  is  always  fulfilled. 

Based  on  an  idea  set  forth  by  Blackstock  [1]  in  the  realm  of 
non-linear  acoustics,  Szabo  proposed  a  way  to  include 
attenuation  and  dispersion  effects  directly  in  the  time  domain 
for  both  non-linear  [2]  and  linear  propagation  [3-5]  in  linear 
media  through  the  inclusion  of  the  so-called  causal 
convolutional  propagation  operator.  By  deriving  a  time 
domain  version  [3]  of  the  Kramer-Kronig  relationships  (K-K), 
[6]  he  arrived  at  a  general  form  for  the  operator.  Szabo's 
operator  was  originally  defined  in  the  context  of  lossy  media 
obeying  a  frequency  power  law  attenuation.  Waters  et  al.  [7] 
showed  that  Szabo's  operator  could  be  used  for  a  broader  class 
of  media  (as  originally  postulated  by  Blackstock),  provided 
the  attenuation  possess  a  Fourier  transform  in  a  distributional 
sense. 

Norton  and  Novarini  [8]  clarified  the  use  of  the  operator  and 
its  relationship  to  the  time  domain  propagation  factor.  The 
time  domain  propagation  factor  is  the  parameter  that  is 
directly  connected  with  the  dispersive  properties  (attenuation 
and  dispersion)  of  the  medium.  The  modified  wave  equation 
that  contains  the  causal  convolutional  propagation  operator  is 
rewritten  incorporating  the  time  domain  propagation  factor 
explicitly.  The  capability  of  the  time  domain  propagation 
factor  to  correctly  incorporate  the  dispersive  traits  of  the 
medium  into  the  linear  wave  equation  was  verified  by  solving 
the  one-dimensional  scalar,  inhomogeneous  wave  equation  in 
a  dispersive  medium  via  a  fmite-difference-time-domain 
(FDTD)  scheme  [8].  It  was  shown  that  for  propagation  in  an 
isotropic  dispersive  (homogeneous)  medium,  attenuation  and 
dispersion  could  correctly  be  taken  into  account  while  staying 
in  the  time  domain.  Therefore,  if  attenuation  versus  frequency 
is  known  (either  from  measurements  or  theoretically)  and 
posses  a  generalized  Fourier  transform,  then  the  time  domain 
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propagation  factor  (and  so  the  corresponding  causal 
convolutional  propagation  operator)  can  be  obtained.  And 
hence  causal  propagation  can  be  modeled  directly  in  the  time 
domain. 

Norton  and  Novarini  [9]  have  extended  the  technique  to  non¬ 
isotropic  media  (i.e.  media  in  which  the  dispersive  effects  vary 
spatially).  Excellent  agreement  between  the  FDTD  with 
dispersion  and  the  analytic  result  were  observed  for  both  the 
forward  and  backscattered  field.  Enhancements  to  the 
numerical  code  included  making  all  difference  approximations 
of  the  derivatives  (both  space  and  time)  fourth  order  accurate 
as  well  as  making  the  absorbing  boundary  conditions  (utilizing 
the  Complementary  Operator  Method)  fourth  order  accurate 
[10].  Norton  and  Novarini  [11]  have  utilized  this  code  to 
investigate  the  scattering  from  and  propagation  though  bubble 
clouds  in  the  ocean.  The  results  show  that  the  presence  of  the 
bubble  clouds  introduces  a  great  deal  of  ringing  that  stays  near 
the  surface.  And  that  there  can  be  a  scattering  effect  from  the 
bubble  cloud  before  the  scattering  event  at  the  air/water 
interface,  resulting  in  returns  from  the  bubble  cloud  arriving  at 
a  receiver  prior  to  the  interface  scattering  return.  Norton  and 
Novarini  [12]  have  also  applied  this  code  to  scattering  and 
propagation  through  biological  tissues,  showing  the  effect  that 
dispersion  has  on  the  backscattered  signal.  The  target 
composition  varied  from  brain,  heart,  kidney,  liver,  and  tendon. 
Norton  [13,14]  has  extended  the  technique  to  include 
heterogeneous  dispersive  media,  recasting  the  modified  wave 
equation  so  that  it  has  a  spatially  varying  bulk  modulus  and 
density.  And  finally  Norton  and  Purrington  [15]  have  replaced 
the  loss  term  found  in  the  Westervelt  equation  with  the  time 
domain  propagation  factor.  They  showed  that  by  employing 
the  time  domain  propagation  factor  the  full  dispersive 
characteristics  of  the  medium  are  properly  incorporated  into 
the  solution,  which  was  lacking  using  the  traditional  loss  term. 

The  paper  is  divided  into  the  following  sections.  First,  the 
linear  wave  equation  for  sound  in  fluids  with  inhomogeneous 
flow  will  be  developed.  Next,  the  introduction  of  dispersive 
effects  into  the  previous  developed  wave  equation  will  be 
accomplished  through  the  development  and  inclusion  of  the 
causal  time  domain  propagation  factor.  Followed  by  a  one¬ 
dimensional  example,  showing  proof  of  concept.  The  resulting 
modified  linear  wave  equation  is  solved  via  the  method  of 
finite  differences. 

II.  Wave  equation  for  sound  in  fluids  with  unsteady 

INHOMOGENEOUS  DISPERSIVE  FLOW 
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The  field  variables  are  expressed  as  sums  of  ambient 
quantities  (Subscript  0)  and  acoustic  perturbations  (primed 
quantities)  so  that  the  linearized  acoustic  equations  can  be 
formed.  That  is 


P=Po  +  P'(x,  0 
P=Po+ff(x>  0 
v=  v0  +  v'(x,f) 
s-  s0  +  s'(x,t). 
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Where,  Dt  =  —  +  V0  •  V  is  the  time  derivative  following 
at 

the  ambient  flow.  Simplification  of  the  above  linearized 
equations  results  when  neglecting  terms  of  2nd  order  and 
higher  in  1/F  and  1/T,  where  F  is  the  length  scale  over  which 
ambient  field  quantities  have  appreciable  spatial  variation  and 
T  is  the  corresponding  time  scale.  This  along  with  the 
introduction  of  a  potential  results  in 


P  =  -pDtO, 


(4) 


A.  Unsteady  Inhomogeneous  Flow 

The  following  development  is  based  on  the  work  by  Pierce 
[16]  and  Godin  [17].  Starting  with  the  full  non-linear  fluid- 
dynamic  equations  for  compressible  fluid  of  uniform 
composition  in  the  absence  of  dissipation  can  be  written  as 
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v'  =  -V®  +  o(l')+0(t')  (5) 

where  is  the  acoustic  part  of  the  flow  velocity 
perturbation.  One  finally  obtains  the  following  expression 
describing  acoustic  propagation  in  non-dispersive  unsteady 
inhomogeneous  flow, 


— V«  (c?VO)-Dt  ^Dt<D  =0 
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where  c  is  a  reference  velocity  and  p  is  the  ambient  density. 

B.  Time  Domain  Propagation  Factor 

The  wave  equation  Eq.  (6)  developed  in  the  last  section 
does  not  account  for  a  dispersive  medium.  Utilizing  the 
approach  as  set  forward  by  Szabo  [3]  an  additional  term  will 
be  introduced  to  the  wave  equation  Eq.  (6)  developed  in  the 
last  section.  The  details  of  the  development  of  the  time 
domain  propagation  factor  can  be  found  in  Ref  [8].  The 
resulting  wave  equation  appropriate  for  a  medium  exhibiting 
dispersive  unsteady  inhomogeneous  flow  is  the  following, 
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I  k(x,z )  dt 


=  S(x-xs)S(z~zs)s(t) 


where  the  spatial  and  temporal  dependence  is  explicitly 
written.  The  bulk  modulus  of  the  medium  is  denoted  by 
and  p{x,z)  is  the  density.  The  two  deltas  on  the 
RHS  are  Dirac  delta  functions  and  s(t)  is  the  source  signature 
at  the  source  location  (xs,zs).  The  time  domain  propagation 
factor  r(r)  is  equal  to 


A.  Simple  Nepheloid  Layer 

In  this  section  a  numerical  experiment  will  be  carried  out  in 
two-dimensions  (X  -  horizontal  axis  and  Z  -  vertical  axis). 
The  environment  consists  of  45m  deep  static  water  layer  over 
lying  a  15m  deep  dispersive  moving  layer  overlying  a  10m 
deep  dispersive  sediment  layer.  Figure  1  depicts  the 
environment.  The  source  consisting  of  ten  periods  of  a  5kHz 
sinusoidal  burst  that  is  Gaussian  weighted  is  located  at  a  range 
of  50m  and  lm  above  the  sediment  layer  (59m  depth).  The 
source  is  located  within  the  dispersive  nepheloid  layer.  Four 
receivers  are  also  located  within  this  layer,  located  at  the 
following  locations,  R1  (10m,  59m);  R2  (25m,  59m);  R3  (75m, 
59m);  R4  (90m,  59m).  Note  that  R1  and  R4  are  equidistance 
from  the  source  and  R2,  R3  likewise  are  equidistance  from  the 
source.  The  reference  sound  speed  in  the  water,  including  the 
dispersive  nepheloid  layer,  is  1500m/s  and  in  the  sediment, 
1700m/s.  The  density  in  the  non-moving  water  layer  is 
1000kg/m3  and  in  the  dispersive  nepheloid  layer  it  varies 
linearly  from  1000kg/m3  to  1050kg/m3.  The  density  in  the 
sediment  is  2000kg/m3. 
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Figure  1.  Graphic  depiction  of  environment. 


r(V)  =  -2  1 +(r)FT-l{a(co)}  (8) 

where  Ot{co)  is  the  attenuation  as  a  function  of  frequency. 
The  attenuation  and  phase  velocity  are  causally  connected 
through  a  Hilbert  transform.  Thus  including  the  attenuation 
will  also  include  the  phase  velocity  in  the  calculation.  The  step 
function  1+(Y)  is  given  as 

0  T  <  0 

l  +  (d  =  j  7-=0  (9) 

1  T>  0. 

It  should  be  noted  that  Szabo’ s  original  use  of  the  operator 
was  limited  to  a  power  law  attenuation  form,  excluding  some 
values  of  the  exponent.  Waters  et  al.  [7],  using  distributional 
analysis,  showed  that  any  attenuation  form  can  be  used  as  long 
as  an  associated  causal  phase  velocity  exists. 

III.  Numerical  Example 


B.  Attenuation 

The  attenuation  is  assumed  to  consist  of  a  power  law 
dependence  on  frequency,  that  is, 

a(f)=kfy ■  (io) 

Where  k  is  a  constant  dependent  upon  the  medium.  For  sand 
it  is  0.3(dB/m/kHz).  And  the  exponent  can  vary  between  1  and 
2.  The  dispersive  nepheloid  layer  is  divided  into  three  equal 
layers  of  5m  thickness  each.  It  is  assumed  that  the  density  of 
resuspended  sediments  increases  linearly  with  depth  thus  the 
dispersive  characteristics  will  likewise  vary  with  depth. 
Additionally  the  resuspended  sediment  will  assume  to  be  silty- 
clay.  Thus  three  different  attenuation  profiles  will  be  modeled, 
using  progressively  larger  values  of  the  constant  k.  The  values 
for  the  constant  are  0.01  (dB/m/kHz),  0.05  (dB/m/kHz),  and 
0.075  (dB/m/kHz).  A  linear  dependence  on  frequency  is 
assumed  for  all  cases,  (i.e.  y=l).  Figure  2  depicts  the 
attenuation  verses  frequency  for  the  dispersive  nepheloid  layer 
and  Figure  3  depicts  the  associated  causal  phase  velocity. 
Figure  4  and  5  depicts  the  attenuation  and  phase  velocity 
verses  frequency  for  the  sediment  layer.  Thus  there  will  be 
three  different  time  domain  propagation  factors  describing  the 
dispersive  characteristics  of  the  dispersive  nepheloid  layer  and 
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one  time  domain  propagation  factor  describing  the  dispersive 
sediment  layer. 
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Figure  5.  Phase  Velocity  vs.  Frequency  for  the  Sediment 
Layer. 


Figure  2.  Attenuation  vs.  Frequency  for  the  Three  Depth 
Segments  of  the  Nepheloid  Layer. 
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Figure  3.  Phase  Velocities  vs.  Frequency  for  the  Three 
Depth  Segments  of  the  Nepheloid  Layer. 


Figure  4.  Attenuation  vs.  Frequency  for  the  Sediment  Layer. 


C.  Flow  Velocities 

The  two  dimensional  flow  velocities  used  within  the 
dispersive  nepheloid  layer  is  based  on  measurements  collected 
during  a  2009  experiment  conducted  in  the  Gulf  of  Mexico 
[18].  Figure  6  and  Figure  7  depicts  the  X  and  Z  component  of 
the  flow  velocity.  The  positive  X-axis  direction  is  to  the  right 
while  the  positive  Z-axis  direction  is  downward.  As  can  be 
seen  from  Figure  6,  the  water  flow  is  to  the  left  and  from 
Figure  7  the  water  will  vary  in  the  vertical  direction  based 
upon  the  depth  within  the  dispersive  nepheloid  layer. 


Figure  6.  X-  Axis  Component  of  Fluid  Velocity 


Figure  7.  Z-  Axis  Component  of  Fluid  Velocity 


IV.  Results 


The  source  is  turned  on  and  the  times  series  is  collected  at 
the  four  receiver  locations.  Figure  8  depicts  the  signal 
amplitude  at  receivers  R1  (solid  black  line)  and  R3  (solid  red 
line).  If  the  time  series  at  R2  and  R4  were  also  plotted,  they 
would  be  superimposed  upon  R1  and  R3  since  the  flow 
velocities  are  so  low  and  the  propagation  distances  are  small. 
However  to  observe  the  difference  in  the  time  series  of  Rl,  R4 
and  R2,  R3  they  were  differenced  and  then  Fourier 
Transformed  to  the  frequency  domain  and  then  plotted.  Figure 
9  depicts  the  differences  between  the  time  series  at  Rl,  R4 
(solid  black  line)  and  R2,  R3  (dotted  black  line).  First  note 
that  the  amplitude  is  different  for  both  sets.  This  is  because  of 
the  different  propagation  distances.  The  distance  is  less  to  R2 
and  R3  thus  the  signal  is  stronger  (i.e.  suffered  less  attenuation) 
than  the  signals  collected  at  Rl  and  R4.  Also  note  that  the 
peak  frequency  is  no  longer  5kHz,  it  is  down  shifted  due  to  the 
linear  frequency  dependence  of  the  attenuation  (higher 
frequencies  are  attenuated  more).  The  only  difference  between 
the  collected  time  series  at  Rl,  R4  and  R2,  R3  is  that  for  Rl 
and  R2  the  signal  propagated  with  the  flow  velocities  and  for 
R3,  R4  it  propagated  against  the  flow  velocities.  To  insure  that 
the  numerical  model  was  not  introducing  a  bias  to  the  result, 
the  flow  velocities  were  zeroed  out  and  the  source  once  again 
turned  on.  The  differences  were  again  taken  between  Rl,  R4 
and  R2,  R3  and  are  shown  on  Figure  9  as  a  red  line  and  a  blue 
line  respectively.  Note  that  both  results  show  no  difference 
between  the  time  series.  This  is  to  be  expected  if  the 
environment  is  invariant  with  regard  to  the  propagation 
direction.  Thus  showing  that  there  is  no  bias  being  introduced 
to  the  solution.  Therefore  the  observed  differences  between 
equidistance  receivers  are  due  solely  to  the  presence  of  the 
flow  velocities. 


Figure  8.  Time  Series  at  Receivers  Rl  and  R3. 


Frequency  (Hz) 

Figure  9.  Difference  of  Times  Series  at  Equidistance 
Receivers  from  the  Source. 


V.  Conclusions 

Direct  Time  Domain  modeling  of  a  simple  two-dimensional 
dispersive  nepheloid  layer  has  been  performed  via  the 
numerical  solution  of  a  modified  linear  wave  equation.  It  was 
observed  that  for  even  small  flow  velocities,  the  associated 
time  series  at  receivers  located  equidistance  from  the  source 
showed  appreciable  differences. 
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